#This is the replication file for Thomann, E. van Engen, N. and L. Tummers. Forthcoming. The necessity of discretion: a behavioral evaluation of bottom-up implementation theory.  Journal of Public Administration Research and Theory.

# This code replicates the analysis of dataset 2, recoding method.

#R code based  on Thomann, Eva and Stefan Wittwer (2016). Performing fuzzy- and crisp set QCA with R: A user-oriented beginner's guide. Version 19.9.2016. URL: http://www.evathomann.com/links/qca-r-manual [date of access: 25.11.2016].


#clear workspace
rm(list = ls())


library(lattice); library(arm); library(xtable); library(foreign); library(psych); library(directlabels); library(betareg); library(VIM); library(base); library(plyr); library(dplyr); library(QCA); library(QCAGUI); library(SetMethods)

setwd("C:/Users/Eva/Dropbox/Arbeit/Publikationen/Papers/Two-factor theory/Material/Data/Dataset 2")


##############PREPARING THE DATASET###########

tfr <- read.csv("tfr_d2.csv", header = TRUE, sep = ";",  dec = ",")
describe(tfr)

#convert -99 into missings
tfr[tfr==-99] <- NA
describe(tfr)

#check missings

colnames(tfr)
describe(tfr)

missings <- as.numeric(complete.cases(tfr))
tabmiss <- table(missings)
prop.table(tabmiss)

#identify  cases with missings on all indicators for a causal factor
##for W
colnames(tfr)

tfr$w1miss <- as.numeric(is.na(tfr$W1))
tabw1miss <- table(tfr$w1miss)
round(prop.table(tabw1miss), digits=3)


tfr$w2miss <- as.numeric(is.na(tfr$W2))
tabw2miss <- table(tfr$w2miss)
round(prop.table(tabw2miss), digits=3)

tfr$w3miss <- as.numeric(is.na(tfr$W3))
tabw3miss <- table(tfr$w3miss)
round(prop.table(tabw3miss), digits=3)

tfr$w4miss <- as.numeric(is.na(tfr$W4))
tabw4miss <- table(tfr$w4miss)
round(prop.table(tabw4miss), digits=3)

tfr$w5miss <- as.numeric(is.na(tfr$W5))
tabw5miss <- table(tfr$w5miss)
round(prop.table(tabw5miss), digits=3)

wmiss <- subset(tfr, select = c("w1miss", "w2miss", "w3miss", "w4miss", "w5miss"))
tfr$wmiss <- rowSums(wmiss)
tfr$wmiss
tfr$wmissout <- as.numeric(tfr$wmiss ==5)
tfr$wmissout
tabwmissout <- table(tfr$wmissout)
prop.table(tabwmissout)


##for SP

tfr$sp1miss <- as.numeric(is.na(tfr$SP1))
tabsp1miss <- table(tfr$sp1miss)
round(prop.table(tabsp1miss), digits=3)


tfr$sp2miss <- as.numeric(is.na(tfr$SP2))
tabsp2miss <- table(tfr$sp2miss)
round(prop.table(tabsp2miss), digits=3)

tfr$sp3miss <- as.numeric(is.na(tfr$SP3))
tabsp3miss <- table(tfr$sp3miss)
round(prop.table(tabsp3miss), digits=3)

tfr$sp4miss <- as.numeric(is.na(tfr$SP4))
tabsp4miss <- table(tfr$sp4miss)
round(prop.table(tabsp4miss), digits=3)

tfr$sp5miss <- as.numeric(is.na(tfr$SP5))
tabsp5miss <- table(tfr$sp5miss)
round(prop.table(tabsp5miss), digits=3)

tfr$sp6miss <- as.numeric(is.na(tfr$SP6))
tabsp6miss <- table(tfr$sp6miss)
round(prop.table(tabsp6miss), digits=3)

spmiss <- subset(tfr, select = c("sp1miss", "sp2miss", "sp3miss", "sp4miss", "sp5miss", "sp6miss"))
tfr$spmiss <- rowSums(spmiss)
tfr$spmiss
tfr$spmissout <- as.numeric(tfr$spmiss==6)
tfr$spmissout
tabspmissout <- table(tfr$spmissout)
prop.table(tabspmissout)

##for TP
tfr$tp1miss <- as.numeric(is.na(tfr$TP1))
tabtp1miss <- table(tfr$tp1miss)
round(prop.table(tabtp1miss), digits=3)


tfr$tp2miss <- as.numeric(is.na(tfr$TP2))
tabtp2miss <- table(tfr$tp2miss)
round(prop.table(tabtp2miss), digits=3)

tfr$tp3miss <- as.numeric(is.na(tfr$TP3))
tabtp3miss <- table(tfr$tp3miss)
round(prop.table(tabtp3miss), digits=3)

tfr$tp4miss <- as.numeric(is.na(tfr$TP4))
tabtp4miss <- table(tfr$tp4miss)
round(prop.table(tabtp4miss), digits=3)

tfr$tp5miss <- as.numeric(is.na(tfr$TP5))
tabtp5miss <- table(tfr$tp5miss)
round(prop.table(tabtp5miss), digits=3)

tfr$tp6miss <- as.numeric(is.na(tfr$TP6))
tabtp6miss <- table(tfr$tp6miss)
round(prop.table(tabtp6miss), digits=3)


tpmiss <- subset(tfr, select = c("tp1miss", "tp2miss", "tp3miss", "tp4miss", "tp5miss", "tp6miss"))
tfr$tpmiss <- rowSums(tpmiss)

tfr$tpmissout <- as.numeric(tfr$tpmiss==6)
tfr$tpmissout
tabtpmissout <- table(tfr$tpmissout)
prop.table(tabtpmissout)

##for OP
tfr$op1miss <- as.numeric(is.na(tfr$OP1))
tabop1miss <- table(tfr$op1miss)
round(prop.table(tabop1miss), digits=3)

tfr$op2miss <- as.numeric(is.na(tfr$OP2))
tabop2miss <- table(tfr$op2miss)
round(prop.table(tabop2miss), digits=3)

tfr$op3miss <- as.numeric(is.na(tfr$OP3))
tabop3miss <- table(tfr$op3miss)
round(prop.table(tabop3miss), digits=3)

tfr$op4miss <- as.numeric(is.na(tfr$OP4))
tabop4miss <- table(tfr$op4miss)
round(prop.table(tabop4miss), digits=3)

tfr$op5miss <- as.numeric(is.na(tfr$OP5))
tabop5miss <- table(tfr$op5miss)
round(prop.table(tabop5miss), digits=3)

tfr$op6miss <- as.numeric(is.na(tfr$OP6))
tabop6miss <- table(tfr$op6miss)
round(prop.table(tabop6miss), digits=3)

opmiss <- subset(tfr, select = c("op1miss", "op2miss", "op3miss", "op4miss", "op5miss", "op6miss"))
tfr$opmiss <- rowSums(opmiss)

tfr$opmissout <- as.numeric(tfr$opmiss==6)
tfr$opmissout
tabopmissout <- table(tfr$opmissout)
prop.table(tabopmissout)

##for SM

tfr$sm1miss <- as.numeric(is.na(tfr$SM1))
tabsm1miss <- table(tfr$sm1miss)
round(prop.table(tabsm1miss), digits=3)

tfr$sm2miss <- as.numeric(is.na(tfr$SM2))
tabsm2miss <- table(tfr$sm2miss)
round(prop.table(tabsm2miss), digits=3)

tfr$sm3miss <- as.numeric(is.na(tfr$SM3))
tabsm3miss <- table(tfr$sm3miss)
round(prop.table(tabsm3miss), digits=3)

tfr$sm4miss <- as.numeric(is.na(tfr$SM4))
tabsm4miss <- table(tfr$sm4miss)
round(prop.table(tabsm4miss), digits=3)
colnames(tfr)


smmiss <- subset(tfr, select = c("sm1miss", "sm2miss", "sm3miss", "sm4miss"))
tfr$smmiss <- rowSums(smmiss)

tfr$smmissout <- as.numeric(tfr$smmiss==12)
tfr$smmissout
tabsmmissout <- table(tfr$smmissout)
prop.table(tabsmmissout)

##for CM

tfr$cm1miss <- as.numeric(is.na(tfr$CM1))
tabcm1miss <- table(tfr$cm1miss)
round(prop.table(tabcm1miss), digits=3)


tfr$cm2miss <- as.numeric(is.na(tfr$CM2))
tabcm2miss <- table(tfr$cm2miss)
round(prop.table(tabcm2miss), digits=3)

tfr$cm3miss <- as.numeric(is.na(tfr$CM3))
tabcm3miss <- table(tfr$cm3miss)
round(prop.table(tabcm3miss), digits=3)

tfr$cm4miss <- as.numeric(is.na(tfr$CM4))
tabcm4miss <- table(tfr$cm4miss)
round(prop.table(tabcm4miss), digits=3)


cmmiss <- subset(tfr, select = c("cm1miss", "cm2miss", "cm3miss", "cm4miss"))
tfr$cmmiss <- rowSums(cmmiss)

tfr$cmmissout <- as.numeric(tfr$cmmiss ==4)
tfr$cmmissout
tabcmmissout <- table(tfr$cmmissout)
prop.table(tabcmmissout)

# calculate dropout rate due to missings

dropout <- subset(tfr, select = c("wmissout", "spmissout", "tpmissout", "opmissout", "smmissout", "cmmissout"))
tfr$dropout <- as.numeric(rowSums(dropout) == 6)
tabdropout <- table(tfr$dropout)
prop.table(tabdropout)

#delete cases with missings on all indicators for 1 causal factor

tfr = tfr[tfr$wmissout == 0,]
tfr = tfr[tfr$spmissout == 0,]
tfr = tfr[tfr$tpmissout == 0,]
tfr = tfr[tfr$opmissout == 0,]
tfr = tfr[tfr$smmissout == 0,]
tfr = tfr[tfr$cmmissout == 0,]
tfr

describe(tfr)

#check missings
missings <- as.numeric(complete.cases(tfr))
tabmiss <- table(missings)
prop.table(tabmiss)

#calculate dropout rate - from 1096 to 1087

do <- round(1-(1087/1096), digits=3)
do


#####CALIBRATION########

#calibrate indicator sets, using indirect calibration method
#the values on the likert scales for powerlessness and meaninglessness need to be reverted because they should mean powerfulness and meaningfulness
##indirect method

###for W
# rename raw variables
tfr <- rename(tfr, w1 = W1)
tfr <- rename(tfr, w2 = W2)
tfr <- rename(tfr, w3 = W3)
tfr <- rename(tfr, w4 = W4)
tfr <- rename(tfr, w5 = W5)
colnames (tfr)


tfr$W1 <- NA
tfr$W1 <- recode(tfr$w1, "5 = 1; 4 = 0.8; 3 = 0.4; 2 = 0.2; 1 = 0; else = copy")
describe(tfr$W1)

tfr$W2 <- NA
tfr$W2 <- recode(tfr$w2, "5 = 1; 4 = 0.8; 3 = 0.4; 2 = 0.2; 1 = 0; else = copy")
describe(tfr$W2)

tfr$W3 <- NA
tfr$W3 <- recode(tfr$w3, "5 = 1; 4 = 0.8; 3 = 0.4; 2 = 0.2; 1 = 0; else = copy")
describe(tfr$W3)

tfr$W4 <- NA
tfr$W4 <- recode(tfr$w4, "5 = 1; 4 = 0.8; 3 = 0.4; 2 = 0.2; 1 = 0; else = copy")
describe(tfr$W4)

tfr$W5 <- NA
tfr$W5 <- recode(tfr$w5, "5 = 1; 4 = 0.8; 3 = 0.4; 2 = 0.2; 1 = 0; else = copy")
describe(tfr$W5)

## for SP
# rename raw variables
tfr <- rename(tfr, sp1 = SP1)
tfr <- rename(tfr, sp2 = SP2)
tfr <- rename(tfr, sp3 = SP3)
tfr <- rename(tfr, sp4 = SP4)
tfr <- rename(tfr, sp5 = SP5)
tfr <- rename(tfr, sp6 = SP6)
colnames (tfr)

tfr$SP1 <- NA
tfr$SP1 <- recode(tfr$sp1, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SP1)


tfr$SP2 <- NA
tfr$SP2 <- recode(tfr$sp2, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SP2)


tfr$SP3 <- NA
tfr$SP3 <- recode(tfr$sp3, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SP3)


tfr$SP4 <- NA
tfr$SP4 <- recode(tfr$sp4, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SP4)


tfr$SP5 <- NA
tfr$SP5 <- recode(tfr$sp5, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SP5)

tfr$SP6 <- NA
tfr$SP6 <- recode(tfr$sp6, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SP6)

## for TP
colnames(tfr)
tfr <- rename(tfr, tp1 = TP1)
tfr <- rename(tfr, tp2 = TP2)
tfr <- rename(tfr, tp3 = TP3)
tfr <- rename(tfr, tp4 = TP4)
tfr <- rename(tfr, tp5 = TP5)
tfr <- rename(tfr, tp6 = TP6)
colnames (tfr)

tfr$TP1 <- NA
tfr$TP1 <- recode(tfr$tp1, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$TP1)


tfr$TP2 <- NA
tfr$TP2 <- recode(tfr$tp2, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$TP2)


tfr$TP3 <- NA
tfr$TP3 <- recode(tfr$tp3, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$TP3)


tfr$TP4 <- NA
tfr$TP4 <- recode(tfr$tp4, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$TP4)


tfr$TP5 <- NA
tfr$TP5 <- recode(tfr$tp5, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$TP5)

tfr$TP6 <- NA
tfr$TP6 <- recode(tfr$tp6, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$TP6)


##for OP

tfr <- rename(tfr, op1 = OP1)
tfr <- rename(tfr, op2 = OP2)
tfr <- rename(tfr, op3 = OP3)
tfr <- rename(tfr, op4 = OP4)
tfr <- rename(tfr, op5 = OP5)
tfr <- rename(tfr, op6 = OP6)
colnames (tfr)

tfr$OP1 <- NA
tfr$OP1 <- recode(tfr$op1, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$OP1)


tfr$OP2 <- NA
tfr$OP2 <- recode(tfr$op2, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$OP2)


tfr$OP3 <- NA
tfr$OP3 <- recode(tfr$op3, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$OP3)


tfr$OP4 <- NA
tfr$OP4 <- recode(tfr$op4, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$OP4)


tfr$OP5 <- NA
tfr$OP5 <- recode(tfr$op5, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$OP5)

tfr$OP6 <- NA
tfr$OP6 <- recode(tfr$op6, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$OP6)


##for SM
colnames(tfr)
tfr <- rename(tfr, sm1 = SM1)
tfr <- rename(tfr, sm2 = SM2)
tfr <- rename(tfr, sm3 = SM3)
tfr <- rename(tfr, sm4 = SM4)

colnames (tfr)

tfr$SM1 <- NA
tfr$SM1 <- recode(tfr$sm1, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM1)


tfr$SM2 <- NA
tfr$SM2 <- recode(tfr$sm2, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM2)


tfr$SM3 <- NA
tfr$SM3 <- recode(tfr$sm3, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM3)


tfr$SM4 <- NA
tfr$SM4 <- recode(tfr$sm4, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$SM4)


##For CM
## for CM
tfr <- rename(tfr, cm1 = CM1)
tfr <- rename(tfr, cm2 = CM2)
tfr <- rename(tfr, cm3 = CM3)
tfr <- rename(tfr, cm4 = CM4)
colnames (tfr)

tfr$CM1 <- NA
tfr$CM1 <- recode(tfr$cm1, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$CM1)


tfr$CM2 <- NA
tfr$CM2 <- recode(tfr$cm2, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$CM2)


tfr$CM3 <- NA
tfr$CM3 <- recode(tfr$cm3, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$CM3)


tfr$CM4 <- NA
tfr$CM4 <- recode(tfr$cm4, "5 = 0; 4 = 0.2; 3 = 0.4; 2 = 0.8; 1 = 1; else = copy")
describe(tfr$CM4)

######AGGREGATION#########

#build aggregated sets using the maximum rule, check for values of 0.5, missings, and for skewedness

# for W

tfr$W <- pmax(tfr$W1, tfr$W2, tfr$W3, tfr$W4, tfr$W5, na.rm=TRUE)
tfr$W

checkw <- as.numeric(tfr$W == 0.5)
sum(checkw, na.rm=TRUE)
describe(tfr$W)
sum(is.na(tfr$W))

wskew <- as.numeric(tfr$W > 0.5)
tabw <- table(wskew)
prop.table(tabw)



# for sp

tfr$SP <- pmax(tfr$SP1, tfr$SP2, tfr$SP3, tfr$SP4, tfr$SP5,tfr$SP6,  na.rm=TRUE)
tfr$SP

checksp <- as.numeric(tfr$SP == 0.5)
sum(checksp, na.rm=TRUE)
describe(tfr$SP)
sum(is.na(tfr$SP))

spskew <- as.numeric(tfr$SP > 0.5)
tabsp <- table(spskew)
prop.table(tabsp)

# for tp

tfr$TP <- pmax(tfr$TP1, tfr$TP2, tfr$TP3, tfr$TP4, tfr$TP5,tfr$TP6,  na.rm=TRUE)
tfr$TP

checktp <- as.numeric(tfr$TP == 0.5)
sum(checktp, na.rm=TRUE)
describe(tfr$TP)
sum(is.na(tfr$TP))

tpskew <- as.numeric(tfr$TP > 0.5)
tabtp <- table(tpskew)
prop.table(tabtp)

# for op

tfr$OP <- pmax(tfr$OP1, tfr$OP2, tfr$OP3, tfr$OP4, tfr$OP5,tfr$OP6,  na.rm=TRUE)
tfr$OP

checkop <- as.numeric(tfr$OP == 0.5)
sum(checkop, na.rm=TRUE)
describe(tfr$OP)
sum(is.na(tfr$OP))

opskew <- as.numeric(tfr$OP > 0.5)
tabop <- table(opskew)
prop.table(tabop)


# for sm

tfr$SM <- pmax(tfr$SM1, tfr$SM2, tfr$SM3, tfr$SM4, na.rm=TRUE)
tfr$SM

checksm <- as.numeric(tfr$SM == 0.5)
sum(checksm, na.rm=TRUE)
describe(tfr$SM)
sum(is.na(tfr$SM))

smskew <- as.numeric(tfr$SM > 0.5)
tabsm <- table(smskew)
prop.table(tabsm)

# for cm


tfr$CM <- pmax(tfr$CM1, tfr$CM2, tfr$CM3, tfr$CM4, na.rm = TRUE)
tfr$CM

checkcm <- as.numeric(tfr$CM == 0.5)
sum(checkcm, na.rm=TRUE)
describe(tfr$CM)
sum(is.na(tfr$CM))

cmskew <- as.numeric(tfr$CM > 0.5)
tabcm <- table(cmskew)
prop.table(tabcm)


#exclude cases with set membership 0.5

tfr = tfr[tfr$W != 0.5,]
tfr = tfr[tfr$OP != 0.5,] 
tfr = tfr[tfr$SP != 0.5,]
tfr = tfr[tfr$TP != 0.5,]
tfr = tfr[tfr$SM != 0.5,]
tfr = tfr[tfr$CM != 0.5,]
tfr

describe(tfr)


tff_d2_c2 <- subset(tfr, select = c("W", "SP", "TP", "OP", "SM", "CM"))
write.csv2(tff_d2_c2, file = "tff_d2_c2.csv")




###################ANALYSIS OF NECESSITY############

tff <- read.csv("tff_d2_c2.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)

#analysis of necessity for W and w  

superSubset(tff, outcome = "W", 
            conditions = "OP, SP, TP, SM, CM", 
            incl.cut = 0.9, cov.cut = 0.6)

superSubset(tff, outcome = "~W", 
            conditions = "OP, SP, TP, SM, CM", 
            incl.cut = 0.9, cov.cut = 0.6)

ttw1


#calculate Z-Test for necessary conditions. I apply a RoN threshold of 0.55, so only for NCs that meet this threshold

nw <- 1-tff$W

# for NC1 TP+cm

NC1W <- compute("TP+cm", data=tff)
obsNC1W <- as.numeric(NC1W >= tff$W)

zNC1W <- round(((((sum(obsNC1W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC1W

pzNC1W <- 2*pnorm(-abs(zNC1W))
pzNC1W

# for NC2 TP+sm

NC2W <- compute("TP+sm", data=tff)
obsNC2W <- as.numeric(NC2W >= tff$W)

zNC2W <- round(((((sum(obsNC2W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC2W

pzNC2W <- 2*pnorm(-abs(zNC2W))
pzNC2W

# for NC3 OP+cm

NC3W <- compute("OP+cm", data=tff)
obsNC3W <- as.numeric(NC3W >= tff$W)

zNC3W <- round(((((sum(obsNC3W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC3W

pzNC3W <- 2*pnorm(-abs(zNC3W))
pzNC3W

# for NC4 tp+SM+cm

NC4W <- compute("tp+SM+cm", data=tff)
obsNC4W <- as.numeric(NC4W >= tff$W)

zNC4W <- round(((((sum(obsNC4W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC4W

pzNC4W <- 2*pnorm(-abs(zNC4W))
pzNC4W

# for NC5 sp+TP+SM

NC5W <- compute("sp+TP+SM", data=tff)
obsNC5W <- as.numeric(NC5W >= tff$W)

zNC5W <- round(((((sum(obsNC5W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC5W

pzNC5W <- 2*pnorm(-abs(zNC5W))
pzNC5W

# for NC6 SP+sm+cm

NC6W <- compute("SP+sm+cm", data=tff)
obsNC6W <- as.numeric(NC6W >= tff$W)

zNC6W <- round(((((sum(obsNC6W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC6W

pzNC6W <- 2*pnorm(-abs(zNC6W))
pzNC6W

# for NC7 SP+sm+CM

NC7W <- compute("SP+sm+CM", data=tff)
obsNC7W <- as.numeric(NC7W >= tff$W)

zNC7W <- round(((((sum(obsNC7W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC7W

pzNC7W <- 2*pnorm(-abs(zNC7W))
pzNC7W

# for NC9 SP+tp+cm

NC9W <- compute("SP+tp+cm", data=tff)
obsNC9W <- as.numeric(NC9W >= tff$W)

zNC9W <- round(((((sum(obsNC9W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC9W

pzNC9W <- 2*pnorm(-abs(zNC9W))
pzNC9W

# for NC10 SP+TP+CM

NC10W <- compute("SP+TP+CM", data=tff)
obsNC10W <- as.numeric(NC10W >= tff$W)

zNC10W <- round(((((sum(obsNC10W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC10W

pzNC10W <- 2*pnorm(-abs(zNC10W))
pzNC10W

# for NC11 SP+TP+SM

NC11W <- compute("SP+TP+SM", data=tff)
obsNC11W <- as.numeric(NC11W >= tff$W)

zNC11W <- round(((((sum(obsNC11W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC11W

pzNC11W <- 2*pnorm(-abs(zNC11W))
pzNC11W

# for NC12 op+TP+CM

NC12W <- compute("op+TP+CM", data=tff)
obsNC12W <- as.numeric(NC12W >= tff$W)

zNC12W <- round(((((sum(obsNC12W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC12W

pzNC12W <- 2*pnorm(-abs(zNC12W))
pzNC12W

# for NC13 op+TP+SM

NC13W <- compute("op+TP+SM", data=tff)
obsNC13W <- as.numeric(NC13W >= tff$W)

zNC13W <- round(((((sum(obsNC13W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC13W

pzNC13W <- 2*pnorm(-abs(zNC13W))
pzNC13W

# for NC14 op+sp+TP

NC14W <- compute("op+sp+TP", data=tff)
obsNC14W <- as.numeric(NC14W >= tff$W)

zNC14W <- round(((((sum(obsNC14W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC14W

pzNC14W <- 2*pnorm(-abs(zNC14W))
pzNC14W

# for NC15 op+SP+cm

NC15W <- compute("op+SP+cm", data=tff)
obsNC15W <- as.numeric(NC15W >= tff$W)

zNC15W <- round(((((sum(obsNC15W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC15W

pzNC15W <- 2*pnorm(-abs(zNC15W))
pzNC15W

# for NC16 op+SP+TP

NC16W <- compute("op+SP+TP", data=tff)
obsNC16W <- as.numeric(NC16W >= tff$W)

zNC16W <- round(((((sum(obsNC16W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC16W

pzNC16W <- 2*pnorm(-abs(zNC16W))
pzNC16W

# for NC17 OP+sm+CM

NC17W <- compute("OP+sm+CM", data=tff)
obsNC17W <- as.numeric(NC17W >= tff$W)

zNC17W <- round(((((sum(obsNC17W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC17W

pzNC17W <- 2*pnorm(-abs(zNC17W))
pzNC17W

# for NC18 OP+tp+sm

NC18W <- compute("OP+tp+sm", data=tff)
obsNC18W <- as.numeric(NC18W >= tff$W)

zNC18W <- round(((((sum(obsNC18W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC18W

pzNC18W <- 2*pnorm(-abs(zNC18W))
pzNC18W

# for NC19 OP+TP+SM

NC19W <- compute("OP+TP+SM", data=tff)
obsNC19W <- as.numeric(NC19W >= tff$W)

zNC19W <- round(((((sum(obsNC19W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC19W

pzNC19W <- 2*pnorm(-abs(zNC19W))
pzNC19W

# for NC20 OP+sp+sm

NC20W <- compute("OP+sp+sm", data=tff)
obsNC20W <- as.numeric(NC20W >= tff$W)

zNC20W <- round(((((sum(obsNC20W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC20W

pzNC20W <- 2*pnorm(-abs(zNC20W))
pzNC20W

# for NC21 OP+sp+TP

NC21W <- compute("OP+sp+TP", data=tff)
obsNC21W <- as.numeric(NC21W >= tff$W)

zNC21W <- round(((((sum(obsNC21W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC21W

pzNC21W <- 2*pnorm(-abs(zNC21W))
pzNC21W

# for NC22 OP+SP+SM

NC22W <- compute("OP+SP+SM", data=tff)
obsNC22W <- as.numeric(NC22W >= tff$W)

zNC22W <- round(((((sum(obsNC22W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC22W

pzNC22W <- 2*pnorm(-abs(zNC22W))
pzNC22W

# for NC24 OP+SP+tp

NC24W <- compute("OP+SP+tp", data=tff)
obsNC24W <- as.numeric(NC24W >= tff$W)

zNC24W <- round(((((sum(obsNC24W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC24W

pzNC24W <- 2*pnorm(-abs(zNC24W))
pzNC24W

# for NC25 OP+SP+TP

NC25W <- compute("OP+SP+TP", data=tff)
obsNC25W <- as.numeric(NC25W >= tff$W)

zNC25W <- round(((((sum(obsNC25W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC25W

pzNC25W <- 2*pnorm(-abs(zNC25W))
pzNC25W

# for NC26 op+sp+SM+cm

NC26W <- compute("op+sp+SM+cm ", data=tff)
obsNC26W <- as.numeric(NC26W >= tff$W)

zNC26W <- round(((((sum(obsNC26W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC26W

pzNC26W <- 2*pnorm(-abs(zNC26W))
pzNC26W

# for NC27 op+SP+tp+sm

NC27W <- compute("op+SP+tp+sm", data=tff)
obsNC27W <- as.numeric(NC27W >= tff$W)

zNC27W <- round(((((sum(obsNC27W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC27W

pzNC27W <- 2*pnorm(-abs(zNC27W))
pzNC27W

# for NC28 OP+sp+tp+SM

NC28W <- compute("OP+sp+tp+SM", data=tff)
obsNC28W <- as.numeric(NC28W >= tff$W)

zNC28W <- round(((((sum(obsNC28W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC28W

pzNC28W <- 2*pnorm(-abs(zNC28W))
pzNC28W

# for NC29 op+SP+tp+SM+CM

NC29W <- compute("op+SP+tp+SM+CM", data=tff)
obsNC29W <- as.numeric(NC29W >= tff$W)

zNC29W <- round(((((sum(obsNC29W == 1))/1087) - 0.8)-(1/(2*1087))) / (sqrt((0.8*(1-0.8))/1087)), digits = 1)
zNC29W

pzNC29W <- 2*pnorm(-abs(zNC29W))
pzNC29W





#Test for skewness of necessary conditions. It is tested whether the proportion of cases in all but one 
#specific area of the XY plot is lower than ten per cent.
##First, visual check with XY Plots; then calculate number of cases outside each area; 
#and calulcate percentage outside area; 
#there are a lot of cases on the diagonal, which are attributed to two area. This is why 
#the percentages don't add up to 100.


##for SP + TP + OP

A1NC25W <- round(sum(as.numeric((NC25W <= tff$W) & (NC25W >= nw)))/(1087/100), digits = 1)

A2NC25W <- round(sum(as.numeric((NC25W >= tff$W)&(NC25W >= nw)))/(1087/100), digits = 1)

A3NC25W <- round(sum(as.numeric((NC25W >= tff$W)&(NC25W <= nw)))/(1087/100), digits=1)

A4NC25W <- round(sum(as.numeric((NC25W <= tff$W)&(NC25W <= nw)))/(1087/100), digits = 1)

sum(A1NC25W, A2NC25W, A3NC25W, A4NC25W)

par(mfrow=c(3, 2))
plot(tff$W, NC25W, pch=18,  xlim=c(0, 1), ylim=c(0,1), col="black",
     xaxs="i", yaxs="i",
     main='Skewness test for SP + TP + OP ',
     xlab=' SP + TP + OP ',
     ylab='W')
abline(0,1, col="black")
abline(1, -1, col="black")
text(x=0.5, y=0.8, labels = A1NC25W)
text(x=0.8, y=0.5, labels = A2NC25W)
text(x=0.5, y=0.2, labels = A3NC25W)
text(x=0.2, y=0.5, labels = A4NC25W)


#################ANALYSIS OF SUFFICIENCY FOR W##################

tff <- read.csv("tff_d2_c2.csv", row.names=1, header = TRUE, sep = ";",  dec = ",")
describe(tff)

nw <- 1-tff$W

#analysis of sufficiency for W
# calculate untenable assumptions

nec <- deMorgan("SP + TP + OP", snames = "SP, TP, OP")
nec

#### W1: frequency threshold 5  ####

ttW1 <- truthTable(tff, outcome="W", 
                        conditions = "SP, TP, OP, SM, CM", 
                         n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                        complete=FALSE, show.cases=FALSE)
ttW1

# the high consistency and PRI values suggest that the outcome is highly overdetermined. Check the truth table for the negated outcome:


ttw1 <- truthTable(tff, outcome="~W", 
                   conditions = "SP, TP, OP, SM, CM", 
                   incl.cut=0.9468, n.cut=5, sort.by="incl, n", decreasing=TRUE, 
                   complete=FALSE, show.cases=FALSE)
ttw1

# this truth table has consistency scores but extremely low PRI scores for all rows. 



